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ABSTRACT. Ice velocities observed in 2005/06 at three GPS stations along the Sermeq Avannarleq 
flowline, West Greenland, are used to characterize an observed annual velocity cycle. We attempt to 
reproduce this annual ice velocity cycle using a 1 -D ice-flow model with longitudinal stresses coupled to 
a 1-D hydrology model that governs an empirical basal sliding rule. Seasonal basal sliding velocity is 
parameterized as a perturbation of prescribed winter sliding velocity that is proportional to the rate of 
change of glacier water storage. The coupled model reproduces the broad features of the annual basal 
sliding cycle observed along this flowline, namely a summer speed-up event followed by a fall slowdown 
event. We also evaluate the hypothesis that the observed annual velocity cycle is due to the annual 
calving cycle at the terminus. We demonstrate that the ice acceleration due to a catastrophic calving 
event takes an order of magnitude longer to reach CU/ETH ('Swiss') Camp (46 km upstream of the 
terminus) than is observed. The seasonal acceleration observed at Swiss Camp is therefore unlikely to be 
the result of velocity perturbations propagated upstream via longitudinal coupling. Instead we interpret 
this velocity cycle to reflect the local history of glacier water balance. 


1. INTRODUCTION 

Ice discharge from marine outlet glaciers is a function of 
deformational and basal sliding velocities. It has been 
suggested that relatively small increases in surface ablation 
may result in disproportionately large increases in ice 
discharge via basal sliding (Zwally and others, 2002; 
Bartholomew and others, 2010). Recent interferometric 
synthetic aperture radar (InSAR) observations have con- 
firmed that an annual velocity cycle is spatially widespread 
in the marginal ice of West Greenland. This annual velocity 
cycle is most likely due to seasonal changes in basal sliding 
velocity (Joughin and others, 2008). Recently, however, 
some studies have speculated that projected increases in 
surface meltwater production will likely result in a net 
decrease in basal sliding velocity due to a transition from 
relatively inefficient to efficient subglacial drainage (Schoof, 
2010; Sundal and others, 2011). This motivates the need to 
quantitatively address the physical relation between glacier 
hydrology and basal sliding velocity. We therefore seek a 
computationally efficient means of reproducing the ob- 
served spatial and temporal patterns of basal sliding so that 
we can ultimately explore the likely response of outlet 
glacier dynamics to climate change scenarios. Other obser- 
vations and models indicate that the ice discharge from 
Greenland's marine-terminating glaciers is highly sensitive 
to calving-front perturbations, which are subsequently 
propagated upstream via longitudinal coupling (Holland 
and others, 2008; Joughin and others, 2008; Nick and 
others, 2009). We also explore this alternative mechanism as 
a cause for the observed annual velocity cycle. 


The Sermeq Avannarleq ablation zone in West Greenland 
exhibits an annual ice velocity cycle similar to that of an 
alpine glacier. The critical features of this cycle are a 
summer speed-up event followed by a fall slowdown event 
(Colgan and others, 201 la). A qualitative comparison of this 
annual ice velocity cycle to a modeled annual glacier water 
storage cycle suggests that enhanced (suppressed) basal 
sliding generally occurs during periods of positive (negative) 
rates of change of glacier water storage (Colgan and others, 
201 la). This notion is consistent with alpine glacier studies 
that suggest that changes in basal sliding velocity are due to 
changes in the rate of change of glacier water storage (dS/df 
or the difference between rates of glacier water input and 
output; Kamb and others, 1994; Anderson and others, 2004; 
Bartholomaus and others, 2008, 2011). In this conceptual 
model, three general basal sliding states exist: (1) when local 
meltwater input exceeds the transmission ability of the 
subglacial hydrologic system (i.e. dS/dt>0 or increasing 
glacier water storage); (2) when the transmission ability of 
the subglacial hydrologic system exceeds the local input of 
meltwater (i.e. dS/df <0 or decreasing glacier water storage); 
and (3) when meltwater input and subglacial transmission 
ability are in approximate equilibrium. In alpine settings, 
peak basal sliding velocity can be expected when dS/df 
reaches a maximum (although there is some evidence that 
peak sliding velocity exhibits a slight phase-lag behind peak 
dS/df values; Bartholomaus and others, 2008). Following this 
maximum, both dS/df and basal sliding velocity decrease, 
and dS/df becomes negative during the later part of the melt 
season when efficient conduits can transmit more water than 
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Fig-1 . The terminal 55 km of the Sermeq Avannarleq flowline overlaid on a panchromatic WorldView-1 image (acquired 15 July 2009), with 
distance from the terminus indicated in km. The GC-Net AWSs are denoted in red. Inset shows the location of Sermeq Avannarleq in West 
Greenland. 


is delivered to the englacial and subglacial system through 
surface melt. 

Our goal is to reproduce the annual ice velocity cycle 
observed in the Sermeq Avannarleq ablation zone with 
coupled ice-flow and hydrology models. In this paper, we 
couple a one-dimensional (1-D) (depth-integrated) ice-flow 
model to a 1-D (depth-integrated) glacier hydrology model 
(Colgan and others, 201 1 a) via a semi-empirical sliding rule. 
Our goal is not to reproduce specific observed intra- or 
interannual variations in ice velocity. Rather we attempt to 
reproduce an annual glaciohydrology cycle in dynamic 
equilibrium that reproduces the critical features of the 
observed cycle and may thus serve as a basis for future work 
investigating the influence of interannual variations in 
surface ablation on annual ice displacement. The 530 km 
Sermeq Avannarleq flowline runs up-glacier from its 
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Fig. 2. Observed lOday mean ice surface velocity (grey lines) and 
cumulative positive degree-days (PDD; red lines) in 2005 and 2006 
at Swiss Camp (SC), JAR1 and JAR2 (where available) versus day of 
year. Black lines denote the bi-Gaussian characterization of the 
annual ice surface velocity cycle at each station (Eqn (1)). 


tidewater terminus (km 0 at 69.37° N, 50.28° W) to the main 
ice divide of the Greenland ice sheet (km 530 at 71 .54° N, 
37.81 °W; Fig. 1). This flowline lies within 2 km of three 
Greenland Climate Network (GC-Net; Steffen and Box, 
2001) automatic weather stations (AWSs): JAR2 (km 14.0 at 
69.42° N, 50.08° W), JAR1 (km 32.5 at 69.50° N, 49.70° W) 
and CU/ETH ('Swiss') Camp (km 46.0 at 69.56° N, 49.34° W) 
(all positions reported for 2008). In a companion paper to this 
study (Colgan and others, 2011a) we suggested that in the 
Sermeq Avannarleq ablation zone: (1) englacial water table 
elevation, which may be taken as a proxy for glacier water 
storage, oscillates around levels that are relatively close to 
flotation throughout the year; and (2) observed periods of 
enhanced (suppressed) basal sliding qualitatively correspond 
to modeled periods of increasing (decreasing) glacier water 
storage. In this paper, we propose a semi-empirical and site- 
specific sliding rule that relates variations in the modeled rate 
of change of glacier water storage to observed variations in 
basal sliding velocity. 


2. METHODS 

2.1. Observed annual ice surface velocity cycle 

We characterize the annual ice surface velocity cycle at 
three locations (JAR2, JAR1 and Swiss Camp) along the 
terminal ~55 km of the Sermeq Avannarleq flowline using 
differential GPS observations of 10 day mean ice surface 
velocity in 2005 and 2006 (Larson and others, 2001; Zwally 
and others, 2002). This characterization provides a represen- 
tative annual ice velocity cycle against which the accuracy 
of the modeled annual ice velocity cycle can be assessed. At 
all sites, these observations reveal that the ice moves at 
winter velocity until the beginning of a summer speed-up 
event in which ice velocities increase above winter velocity 
(Fig. 2). The summer speed-up event is followed by a fall 
slowdown event in which ice velocities decrease below 
winter velocities. Using the positive degree-days (PDDs) 
observed at each station as a proxy for melt intensity 
(Ohmura, 2001; Steffen and Box, 2001), the onset of the 
speed-up approximately coincides with the onset of summer 
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melt and reaches a maximum approximately halfway 
through the melt season, while the slowdown occurs after 
the cessation of melt. 

We approximate the annual ice surface velocity cycle 
using two Gaussian curves. These are overlaid on mean 
winter velocity, with one (positive) curve representing the 
summer speed-up event and the other (negative) curve 
representing the fall slowdown event. The use of a bi- 
Gaussian function allows the amplitude, width and timing of 
both the summer and fall events to be parameterized 
independently. Thus, we characterize surface velocity, u s , 
as a function of day of year, j, according to: 


u s = u w + (u max - u w ) exp 


- [u w - (u w - Umin)] exp 


J~h 


J J min 
^min 


0 ) 


Table 1. The value of each parameter in the bi-Gaussian 
characterization of the annual ice surface velocity cycle at JAR2, 
JAR1 and Swiss Camp (Eqn (1)) 



JAR2 

JAR1 

Swiss Camp 

u w (m a ’) 

105 

66 

113 

Umax (m a _l ) 

195 

95 

175 

Umin (m a -1 ) 

80 

59 

101 

jmax (days) 

175 

200 

200 

fmin (days) 

255 

250 

235 

dmax (days) 

30 

25 

12 

cLn (days) 

30 

20 

25 


water storage calculated by the hydrology model into 
variations in basal sliding velocity. 


where u w is the mean winter velocity and u max and u min are 
the summer maximum and fall minimum velocities, respect- 
ively. The remaining four parameters govern the timing and 
shape of the speed-up and slowdown curves: y max (y min ) is the 
day of year of summer maximum (fall minimum) velocity, 
while c/ max (d m \ n ) is the duration of the summer (fall) velocity 
anomaly. All parameters were prescribed manually by visual 
inspection (Table 1). For y' max and y min we assess an estimated 
uncertainty equivalent to the temporal resolution of the 
velocity data (±10 days). This bi-Gaussian characterization 
was fitted to the aggregated 2005 and 2006 velocity data at 
each station (Fig. 2). 

2.2. Ice-flow model 

We apply a longitudinally coupled 1-D (depth- integrated) 
ice-flow model to the Sermeq Avannarleq flowline. This 
model solves for the rate of change in ice thickness (dH/dt) 
according to mass conservation: 


where b is the annual mass balance, Q is the ice discharge 
per unit width and dQ/dx is the horizontal divergence of ice 
discharge. To generate dynamic equilibrium ice geometry 
and velocity fields, the ice-flow model was subjected to a 
1000 year spin-up that was initialized with present-day ice 
geometry and a 'cooler' climate with no hydrology cycle 
(described in Section 2.4). We characterize dynamic equi- 
librium as the transient solution of Eqn (2) that exhibits no 
significant changes in ice thickness (\dH/dt\ < 1 m a -1 ) during 
the last 100 years of spin-up. Alternative approaches would 
be to produce: (1) a fully transient non-equilibrium present- 
day snapshot of flowline ice geometry and velocity by spin- 
up under a prescribed climate scenario; or (2) a steady-state 
solution of flowline ice geometry and velocity under 
imposed spin-up conditions. The former would certainly 
be desirable for modeling future flowline evolution, but is 
sensitive to uncertainties in the prescribed climate forcing. 
The latter requires the implementation of boundary condi- 
tions at both ends of the flowline, and steady-state calving 
flux is not precisely known in this instance due to 
uncertainty in flowline delineation. Following the 1 000 year 
spin-up, an annual basal sliding cycle is introduced via the 
coupled 1-D (depth-integrated) hydrology model. We use a 
semi-empirical three-phase sliding rule (described in Section 
2.3) to convert variations in the rate of change of glacier 


2.2.1. Annual balance 

The annual mass balance of a given ice column is the sum of 
the annual surface accumulation, c s , surface ablation, a s , 
basal accumulation, c b, and basal ablation, at,: 

6 = c s + Fa s + q, + at, (3) 

where F is the hydrologic system entry fraction based on the 
ratio of annual surface accumulation to annual surface 
ablation (Pfeffer and others, 1991; Colgan and others, 
2011a). As F is the fraction of ablation assumed to enter 
the glacier hydrology system and eventually flow out of the 
ice sheet, the quantity 1 - F is the fraction of ablation that 
refreezes and does not leave the ice sheet. This assumes that 
purely supraglacial transport to the margin is negligible; at 
Sermeq Avannarleq, all runoff is expected to drain into the 
englacial system via either crevasses or moulins (McGrath 
and others, 2011). Annual surface accumulation is pre- 
scribed as the observed mean annual value over the period 
1991-2000 (Burgess and others, 2010). Annual accumu- 
lation increases from ~0.25 m at the Sermeq Avannarleq 
terminus to a maximum (~0.5 m) at ~1 00 km upstream and 
decreases again to ~0.25 m at the main flow divide (530 km 
upstream). In the ablation zone, annual surface ablation, a s , 
is taken to be a function of elevation, based on previous 
observations: 

a s 7(^s — -^ela) — ^ela (4) 

where 7 is the present-day ablation gradient with elevation 
(Aa s /Az s ; taken as 0.00372; Fausto and others, 2009), z s is 
the ice surface elevation, z e | a is the equilibrium-line altitude 
(ELA) and a e i a is the annual surface ablation at the ELA (taken 
as 0.4 m). The observed regional ELA was ~1 1 25 m over the 
period 1995-99 (Steffen and Box, 2001) and ~1250m over 
the period 1996-2006 (Fausto and others, 2009). We 
prescribe an ELA of 1125 m, as it is more likely to be 
consistent with the steady-state surface mass-balance forcing 
prior to the highly transient post-1990 period. Annual 
surface abalation is distributed over an annual cycle to 
yield surface ablation rate a s using a sine function to 
represent the melt season solar insolation history (cf. fig. 4 in 
Colgan and others, 2011a). 

In the ice-flow model, we assume that annual basal 
accumulation is negligible (Cb~ 0 m a -1 ) and that submarine 
basal ablation is only significant beneath the floating 
tidewater tongue (e.g. Rignot and others, 2010). We use 
the relative magnitudes of ice and englacial water pressures, 
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Fig. 3. The flow-law parameter values, A calculated according to 
Eqn (9) using the Wisconsin enhancement factor values, E (inferred 
by observed ice thickness), and modeled ice temperature values, 7) 
along the Sermeq Avannarleq flowline. Recommended values for 
E=1 ice are shown for comparison (Paterson, 1994). 


Pi and P w , respectively, to determine which flowline nodes 
are grounded (Pi > P w ) or floating (Pi<P w ) in a given time- 
step. During spin-up and dynamic equilibrium we prescribe 
a constant submarine basal ablation rate of ab = 10ma"' to 
all floating nodes. This prescribed rate is less than the 
contemporary submarine basal ablation rate at Sermeq 
Avannarleq, which is estimated to exceed 25 m a -1 (Rignot 
and others, 2010). The ice-flow model, however, does not 
reproduce a floating tongue when contemporary submarine 
basal ablation rates are imposed for the duration of spin-up. 
As our intent is to reproduce the dynamic equilibrium ice 
geometry that precedes the current rapid transient state, we 
depart from the present-day submarine basal ablation rate. 


2.2.2. Ice discharge 

We include depth-averaged longitudinal coupling stress, 
r xx , as a perturbation to the gravitational driving stress 
derived from the shallow-ice approximation (Van der Veen, 
1987; Marshall and others, 2005): 


r)y f) . 

r=-p ig H^ + 2-(H, xx ) (5) 

where r is total driving stress. Depth-averaged longitudinal 
coupling stress (r xx ) is calculated following the approach 
outlined by Van der Veen (1987). This formulation derives 
longitudinal coupling stress by solving a cubic equation 
describing equilibrium forces independently at each node, 
based on ice geometry and prescribed basal sliding velocity, 
u b : 



+ • ' • r. 


dz s dH 3 H d 2 z s 
dx dx 2 dx 2 



fl dH 1 c>z s \ 1 <9u b 
\5 dx 4 dx ) 2A dx 



( 6 ) 

The depth-integrated longitudinally coupled ice velocity 
due to deformation, 17^, may be derived from the equation 
for horizontal shear rate, dujdz (e.g. Marshall and 


others 2005): 


dud = 
dz 

27\(pi g(z s - 




dz X-i 
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-Pig(z s - 


dx 
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d 

dx 



where dzjdx is ice surface slope along the flowline and n is 
an exponent of 3 in the empirical relation between stress 
and strain rate describing ice rheology (Glen, 1955). 
Integrating Eqn (7) twice in the vertical and dividing by H 
yields 


Ud = 


2 A 
n + 2 


u dZi 

KgH or 




H. 


( 8 ) 


We calculate the flow-law parameter, A, as a function of 
both ice temperature, 7) and thickness, H, using an 
Arrhenius-type relation (Huybrechts and others, 1991): 

A( T, H) = E(H)A 0 ( T) exp (9) 


where A 0 is a coefficient that depends on ice temperature 
(taken as 5.47 x 10 10 Pa -3 a -1 when 7~>263.15 K and 
1 .14 x 1 0 -5 Pa -3 a -1 when 7“<263.15K), Q e is the creep 
activation energy of ice (taken as 139kJmoD 1 when 
7“> 263.1 5 K, and eokjmol" 1 when 7~<263.15K), R is the 
ideal gas constant (8.314J mol -1 K -1 ) and T is the ice 
temperature. At each flowline node, we use the steady-state 
ice temperature at 90% depth, derived from independent 
thermodynamic modeling of the flowline (personal commu- 
nication from T. Phillips and others, 2011), to calculate the 
flow-law parameter. The majority of shear occurs at or below 
this depth. Thus, along-flowline variations in basal ice 
temperature result in along-flowline variations in the flow- 
law parameter. 

We enhance the flow-law parameter by a factor £, to 
account for increased deformation due to the presence of 
relatively soft Wisconsin basal ice. This enhancement factor 
linearly transitions from its prescribed value where ice 
thickness exceeds 650 m, to 1 (i.e. no enhancement) where 
ice thickness is <550 m. At flowline nodes where ice is 
>650 m thick, Wisconsin ice is expected to comprise a 
significant portion of the basal ice (Huybrechts, 1994). We 
assume that the uncertainty associated with the calculated 
values of A is small in comparison to uncertainty associated 
with the Wisconsin flow enhancement factor. We evaluate 
dynamic equilibrium ice geometry and velocity fields 
following spin-up with E ranging between 2 and 4 (Reeh, 
1985; Paterson, 1991), as in situ borehole deformation 
measurements beneath nearby Jakobshavn Isbrae indicate 
E> 1 (Luthi and others, 2002). Under the E= 3 scenario, 
calculated flow-law parameter values range between 
8.2 x10 -18 and 4.3 x 1 0 -16 Pa -3 a -1 . For comparison, the 
recommended unenhanced (E=1) flow-law parameter at 
273 K is 2.1 x 1 0' 16 Pa -3 a -1 (Paterson, 1994; Fig. 3). 

Local ice discharge, Q, is obtained by multiplying the 
sum of basal sliding and depth-averaged deformational 
velocities by ice thickness: 

Q = (u b + u d )H. (10) 


Basal sliding velocity is prescribed via a semi-empirical 
sliding rule described in Section 2.3. 

Following the approach taken in previous ice-sheet 
flowline models (e.g. Van der Veen, 1987; Parizek and 
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Alley, 2004), we neglect lateral velocities and stresses 
stemming from divergence and convergence. While this 
assumption is likely valid in interior regions of ice sheets, it 
is less valid near the ice-sheet margin, where substantial 
divergence and convergence can occur. We acknowledge in 
Section 3 that the failure to account for possible lateral 
effects is a potential factor in the systematic overestimation 
of ice velocities at JAR1 and underestimation of ice 
velocities at the terminus. Another inherent shortcoming of 
any 1-D flowline model is the prescription of the ice-flow 
trajectory. For example, if the flowline length has been 
underestimated in the accumulation zone, modeled ice 
discharge across the equilibrium line will be underestimated 
and the model will require a decrease in surface ablation 
rate to maintain the observed ice geometry. In addition, the 
ice-flow trajectory is unlikely to have been constant through 
time. Striations on the ice surface mapped by Thomsen and 
others (1988) suggest that the flowline through JAR1 station 
likely terminated on land just north of Sermeq Avannarleq in 
1985, rather than flowing past JAR2 to the tidewater 
terminus as shown in Figure 1. The recent acceleration of 
Jakobshavn Isbrae, immediately south of Sermeq Avannarleq, 
has caused substantial reorientation of ice flow throughout 
the Sermeq Avannarleq ablation zone (Colgan and others, 
201 1 b). As the flowline used in this study was derived from a 
2005/06 InSAR ice surface velocity field (Joughin and others, 
2010), obtained after the onset of the reorganization of ice 
flow (~1 997; Colgan and others, 2011b), it does not 
accurately reflect the long-term ablation zone flowline 
trajectory. 

2.3. Three-phase basal sliding rule 

The sliding rules employed in glacier models have improved 
with advances in the conceptualization of basal sliding. 
Initial sliding rules prescribed basal sliding velocity as 
proportional to driving stress on the assumption that higher 
driving stress results in greater till deformation (Weertman, 
1957; Kamb, 1970). Observations that subglacial water was 
capable of enhancing ice velocities by lubricating and 
pressurizing the subglacial environment led to parameter- 
izations in which basal sliding velocity was taken as 
proportional to subglacial water pressure (Iken, 1981; Iken 
and others, 1 983), as well as sliding rules that included both 
effective water pressure (P,-P w ) and driving stress (Bind- 
schadler, 1983; Iken and Bindschadler, 1986). Recent 
models have utilized basal sliding rules that are Coulomb 
friction analogues, whereby basal drag is parameterized to 
take on a maximum value at low sliding velocities and 
decrease with decreasing effective pressure and increasing 
sliding velocity (Schoof, 2005; Gagliardini and others, 2007; 
Pimentel and Flowers, 2010). To have predictive value, a 
basal sliding rule should ideally be capable of reproducing 
observed sliding velocities from first principles of hydrology 
and force balance with a minimum of free parameters. FHere 
we describe an empirical and site-specific three-phase basal 
sliding rule that focuses on the hydrologic aspect of basal 
sliding. We regulate the magnitude and sign of a perturb- 
ation to background basal sliding velocities using modeled 
rates of change of glacier water storage. 

The Swiss Camp GPS data indicate that ice velocity is on 
average ~14% faster during the winter than in the midst of 
the slowest motion of the year, during the fall slowdown 
(113 and 99 m a -1 , respectively; Fig. 2; Table 1). This fall 
velocity minimum matches very well the ice surface velocity 
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Fig. 4. Background basal sliding velocity, Ubo, estimated at Swiss 
Camp (SC), JAR1 and JAR2 as the difference between mean winter 
(u w ) and fall minimum (u min ) velocities (Table 1). A least-squares 
linear interpolation/extrapolation to the terminal 65 km of the 
flowline is also shown. 


predicted by internal deformation alone. According to the 
shallow-ice approximation, the ice surface velocity solely 
due to deformation may be calculated as (Hooke, 2005) 

*■<*> = m («»)"""' <"> 

Taking ice thickness, H, as 950 m, ice surface slope, dzjdx, 
as 0.01 and A as 3.33 x 1 0 -16 Pa -3 a -1 (the flow-law 
parameter for ice at 272 K, which is the pressure-melting 
point beneath ~1 km of ice, enhanced by a Wisconsin factor 
of 3), yields an ice surface velocity of 99 m a -1 . 

We interpret this as suggesting that Swiss Camp 
experiences significant background basal sliding velocity 
during the winter, which is suppressed during the fall 
velocity minimum. Thus, at Swiss Camp (and similarly at 
JAR1 and JAR2), the background basal sliding velocity may 
be approximated as the difference between observed mean 
winter (u w ) and fall minimum (u mm ) velocities. We linearly 
interpolate these basal sliding velocities along the flowline 
to provide the background basal sliding boundary condi- 
tion used in spin-up of the ice-flow model (Fig. 4). The 
year-round persistence of the englacial hydrologic system 
in the Sermeq Avannarleq ablation zone provides a 
mechanism capable of maintaining favorable basal sliding 
conditions year-round (i.e. available liquid water; Catania 
and Neumann, 2010). After spin-up, we overlay an annual 
basal sliding velocity cycle on the background basal 
sliding velocity. 

Following theoretical developments in alpine glacio- 
hydrology (Kamb and others, 1994; Anderson and others, 
2004; Bartholomaus and others, 2008, 201 1), we propose a 
sliding rule that depends on the sign of the rate of change of 
glacier water storage to prescribe 'speed-up' during periods 
of increasing glacier water storage and 'slowdown' during 
periods of decreasing glacier water storage. We take rate of 
change of englacial water table elevation (or head, dhjdt) 
as a surrogate for rate of change of glacier water storage 
(Colgan and others, 2011a). We formulate a three-phase 
basal sliding rule that imposes: (1) background basal sliding 
velocity during the winter when dh e /dt&0; (2) enhanced 
basal sliding during positive rates of change of glacier water 
storage {dh e /dt>0); and (3) suppressed basal sliding during 
negative rates of change of glacier water storage (dhjdt< 0). 
We accomplish this by conceptualizing u b as the sum of 
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Fig. 5. Along-flowline distributions of the parameters used to 
calculate basal sliding perturbation (Eqn (13)): (a) sliding rule 
coefficient, k ; (b) subglacial conduits nT 1 in the across-flow 
direction, n c (Colgan and others, 2011a); and (c) daily profiles of 
rate of change of glacier water storage, dhjdt, over an annual cycle 
(Colgan and others, 2011a). Vertical dashed lines denote the 
locations of JAR2, JAR1 and Swiss Camp (SC). 


background basal sliding velocity, u^o, and a perturbation, 
A u h : 

Ub = Ubo + Aub limit: Ub > 0 (12) 

The alpine glaciohydrology literature suggests that we may 
expect basal sliding velocity to scale nonlinearly with the 
rate of change of glacier water storage (i.e. Aub<x(dh e /dt) m , 
where m> 1; Anderson and others, 2004; Bartholomaus and 
others, 2008). We impose m= 3 and express the perturb- 
ation to the background basal velocity as 

Aub = { (nc{x)k{x)\- t|) m sign (t) if m >0.25 m d> 

\o if |^| <0.25 md"' 

(13) 

where n c (x) is the number of subglacial conduits rrT 1 in the 
across- ice-flow direction (Colgan and others, 2011a) and 
k(x) is a tunable site-specific sliding coefficient. We impose 
an arbitrary threshold of 0.25 md -1 (constant in time and 
space) that rates of change of glacier water storage must 
exceed in order to either enhance or suppress background 
basal sliding. We tune the sliding coefficient values (which 
range between ~0.25 and 0.75 m 1/3 d 2/3 ) to reach a 
minimum in the vicinity of JAR1, based on the observation 
that the annual velocity cycle at JAR1 is damped in 
comparison to JAR2 and Swiss Camp (Fig. 5). 

The number of conduits per meter in the across-ice-flow 
direction reflects variations in the configuration of the 
subglacial drainage system with distance upstream, from 
relatively large widely spaced conduits near the terminus to 
relatively small closely spaced conduits near the equilibrium 
line. These differences in subglacial hydrologic system 
configuration can be expected to result in differing sliding 
responses to a given rate of change of glacier water storage. 
For example, the enhanced sliding due to a given meltwater 
input is expected to be greater for a subglacial network 
comprising 'cavities' than for a subglacial network compris- 
ing 'channels' (Schoof, 2010). As the sign of the rate of 
change of glacier water storage, dhjdt, changes between 
positive and negative, it effectively modulates the sign of Aub 


(i.e. specifying whether the perturbation is acting to enhance 
or suppress background basal sliding velocity). By para- 
meterizing Aub so that it goes to zero when \dh e /dt\ falls 
below a critical threshold (taken as 0.25 md -1 ), Eqn (13) 
provides the framework for three phases of basal sliding. 

2.4. Input datasets and boundary conditions 

Following Colgan and others (2011a), the ice-flow model is 
initialized with observed ice surface and bedrock topog- 
raphy (Bamberand others, 2001; Scambos and FHaran, 2002; 
Plummer and others, 2008). During spin-up, the observed 
surface mass balance was adjusted by decreasing surface 
ablation in order to reproduce the observed present-day ice 
geometry. We evaluate dynamic equilibrium ice geometry 
and velocity fields following spin-ups by decreasing surface 
ablation by a factor of 25-75%. Previous Greenland ice 
sheet modeling studies have implemented similar surface 
mass-balance corrections during spin-up in order to achieve 
equilibrium present-day ice geometries (i.e. Fluybrechts, 
1994; Ritz and others, 1997; Parizek and Alley, 2004). This 
adjustment is typically justified by the notion that the 
present-day Greenland ice sheet geometry reflects colder 
'glacial' conditions (i.e. less surface ablation with no change 
in accumulation). Both the observed surface accumulation 
(Burgess and others, 2010) and ablation (Fausto and others, 
2009) datasets have been validated by in situ observations, 
including the GC-Net AWSs along the Sermeq Avannarleq 
flowline (Steffen and Box, 2001; Fig. 1). 

The differential equations describing the evolution of ice 
thickness (Eqns (2-10)) were discretized in space using first- 
order finite volume methods with grid spacing Ax=500m. 
The semi-discrete set of coupled ordinary differential equa- 
tions at the computational nodes was then solved using 
'ode15s', the stiff differential equation solver in MATLAB 
R2008b. During the 1000 year spin-up, the ice-flow model 
was solved with a time-step, At, of 10 years. Following spin- 
up, the ice-flow model was solved concurrently with the 
hydrology model with a 2 day time-step. Post-spin-up, basal 
sliding velocity was calculated according to the three-phase 
sliding rule described above. We apply a second-type (zero 
flux) Neumann boundary condition at the main ice-flow 
divide as the upstream boundary condition (dzjdx= 0 and 
Q= 0 at x= 530 km). When the glacier terminus reaches the 
downstream boundary of the model domain (kmO), the 
terminal node ice discharge (Q te rm or calving flux) is 
calculated as the difference between the ice discharge of 
the adjacent upstream node and the annual mass balance 
(negative in the ablation zone) of the terminal node 
(Qterm = Qterm-i + £>Ax at x=0km). When the glacier ter- 
minus does not reach the downstream boundary of the 
model domain, no calving flux is imposed, as a dynamic 
equilibrium has been achieved in which ice inflow across 
the grounding line is balanced by total ablation (both surface 
and submarine) in the floating ice tongue. 

The 1-D (depth-integrated) hydrology model tracks 
glacier water storage and discharge through time. Glacier 
water input is prescribed based on observed ablation rates, 
whereas glacier water output occurs through conduit 
discharge. Conduit discharge varies in response to the 
dynamic evolution of conduit radius. When coupled, the 
ice-flow and hydrology models receive the same surface 
ablation forcing at each time-step. The ice-flow model 
updates ice geometry used by the hydrology model each 
time-step, while the hydrology model updates the subglacial 
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head values used by the basal sliding rule. Specific 
parameterizations used in the hydrology model are shown 
in Table 2 (cf. animation 1 in Colgan and others, 201 la). 


3. RESULTS 

The GPS-observed ice velocities provide insight into the 
relative importance of basal sliding in the Sermeq Avannar- 
leq ablation zone (Fig. 2). The net annual displacements 
observed at JAR2, JAR1 and Swiss Camp are about 114, 69 
and 115 m, respectively. The fall minimum velocities, taken 
to represent purely deformational velocities (i.e. no basal 
sliding), observed at each of these stations are about 80, 59 
and 101 ma" 1 , respectively. Therefore, basal sliding appears 
to be responsible for about 30, 1 4 and 1 2% of the annual net 
displacement, respectively. These observations suggest that 
summer speed-up events are responsible for about 1 3, 4 and 
4 m of net displacement at JAR2, JAR1 and Swiss Camp, 
respectively, equivalent to about 11,6 and 3% of the annual 
displacement. These increases in net displacement due to 
summer speed-up events are partially offset by decreases in 
net displacement due to fall slowdown events. The fall 
slowdown events suppress annual displacement by about 4, 
1 and 2 m at JAR2, JAR1 and Swiss Camp, respectively, 
equivalent to about 4, 1 and 2% of the annual displacement. 
Thus, at Sermeq Avannarleq, year-round (or 'background') 
basal sliding appears to contribute a larger fraction of annual 
net displacement than seasonal basal sliding. 

We compared modeled dynamic equilibrium ice geom- 
etry and velocity fields of the ice-flow model with observed 
ice surface elevation (Scambos and Haran, 2002) and 
velocity (Joughin and others, 2008) following a 1000 year 
spin-up. We explored a range of Wisconsin enhancement 
factors (2-4) and surface ablation perturbations (decreases of 
25-75% below contemporary rates). While this parameter 
space produces a relatively narrow range of ice surface 
geometries that are similar to the observed profile, it 
produces a relatively wide range of ice velocities that can 
differ substantially from the observed profile (Fig. 6). A 


Table 2. Specific parameterization of the 1-D (depth-integrated) 
hydrology model (notation follows Colgan and others, 201 la) 


Variable 

Definition 

Value 

a 

Glacier hydrology length scale 

20 km 

<t> 

Bulk ice porosity 

0.01 

F r 

Firn meltwater retention fraction 

0.5 

Qg 

Geothermal flux 

57 mW m -2 

f 

Conduit friction factor 

0.05 

pterm 

Conduit spacing at terminus 

0.005 m -1 

..term 

'max 

Maximum conduit radius at terminus 

2m 


comparison between modeled mean ice surface elevation 
and velocity versus observed mean ice surface elevation 
(688 m) and velocity (1 1 6 m a -1 ) along the terminal 50 km of 
the flowline suggests that observations are most accurately 
reproduced with a 50% decrease in surface ablation below 
contemporary rates and a Wisconsin enhancement factor of 
E= 3 (Table 3). This depression in contemporary surface 
ablation may represent a combination of: (1) an under- 
estimation of historical accumulation rate; (2) the increase in 
surface ablation that has occurred since the termination of 
the last glaciation; and (3) error in our delineation of 
the flowline. 

Under all spin-up scenarios, the ice-flow model over- 
estimates ice velocities in the km 25-35 portion of the 
flowline and underestimates ice velocities in the terminal 
~5 km. We speculate that both systematic errors are artifacts 
of the 1-D (flowline) character of the ice-flow model. A 1-D 
model is expected to overestimate ice discharge (and hence 
ice velocity) in reaches of divergent ice flow. Thus, the 
overestimation of ice velocities in the vicinity of JAR1 may 
suggest that divergent ice flow is occurring in this region. 
Similarly, a 1-D model is expected to underestimate ice 
velocities in regions of convergent ice flow. Failure to 
account for convergent flow at the terminus has likely 
contributed to the severe underestimation of terminus ice 
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Fig. 6. (a) Modeled ice surface elevation, z s , and (b) modeled ice surface velocity, u s , over a range of Wisconsin enhancement factors, E, and 
fractional contemporary surface ablation values, a s . Observed ice surface elevation (Scambos and Haran, 2002) and velocity (Joughin and 
others, 2008) are shown for comparison. The basal boundary conditions (BC) for ice surface elevation and velocity are observed bedrock 
elevation and prescribed background basal sliding velocity. Vertical dashed lines denote the locations of JAR2, JAR1 and Swiss Camp (SC). 
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Model time: 363.8 



Distance upstream (km) 

Animation 1 . Animation of the annual glaciohydrology cycle, (a) Surface ablation rate, a s . (b) Bedrock elevation (brown), transient ice 
geometry (black line) and transient englacial water table elevation (blue) (Colgan and others, 201 1 a), (c) Basal sliding velocity calculated 
from the semi-empirical three-phase sliding rule (Eqns (12) and (13)). (d) Bedrock elevation (brown) and transient ice geometry (black 
line) with contour shading to denote ice surface velocity, u s (color bar saturates at 200 ma -1 ). Vertical dashed lines denote the locations 
of JAR2, JAR1 and Swiss Camp (SC). Model time is given in day of year. 

Full movie available at htt p://www.igsoc.org/hyperlink/11l081 A ni mation l .mov. 


thickness. Observations suggest the bedrock overdeepening 
at the terminus contains grounded ice ~600m thick 
(Plummer and others, 2008) while the modeled terminus is 
floating with ice thicknesses of ~300m. 

Generally, however, the ice-flow model achieves a 
dynamic equilibrium ice geometry and velocity field that 
closely matches observations throughout the Sermeq Avan- 
narleq ablation zone (i.e. terminal 60 km of the flowline). 
This suggests that the basal sliding boundary condition 
imposed during the 1000 year spin-up (i.e. 'background 
basal sliding'; Fig. 4) captures the essence of the basal 
sliding profile beneath the terminal ~60 km of the flowline. 
When we allow the annual hydrologic cycle to operate, 
producing seasonal variations in basal sliding around this 
background basal sliding profile, the three-phase basal 
sliding rule produces a reasonable annual basal sliding 
velocity cycle throughout the Sermeq Avannarleq ablation 
zone (Animation 1). While the absolute magnitude and 
seasonal timing of modeled maximum and minimum 
velocities do not precisely match observations (discussed 


in Section 4.1 ), the coupled model produces an annual basal 
sliding cycle that captures the essence of the observed 
annual velocity cycle, namely: (1) prescribed background 
basal sliding velocity during the winter, (2) a summer speed- 
up event, followed by (3) a fall slowdown event. 

4. DISCUSSION 

We developed a representation of basal sliding velocity that 
depends on the local rate of change of glacier water storage. 
Incorporating this basal boundary condition in a 1-D ice- 
flow model captures the essence of the annual velocity cycle 
of the Sermeq Avannarleq flowline (i.e. summer speed-up 
and fall slowdown events). The mean annual glaciohydrol- 
ogy cycle reproduced by this coupled model may serve as a 
basis for future investigations into the transient response of 
the Sermeq Avannarleq flowline to predicted increases in 
surface meltwater production. While we do not drive our 
basal sliding rule with absolute head values, h e/ but rather 
with the rate of change of head, dhjdt, we achieve a 


Table 3. Discrepancy between modeled and observed mean ice surface elevation, Sz s , and mean ice surface velocity, Su s , along the terminal 
50 km of the flowline, under various scenarios of fraction of contemporary surface ablation rate, a s , and Wisconsin flow-law enhancement 
factor, E. Discrepancies are expressed in both absolute and relative values 


E 

<5z s 

m 

0.25 a s 

6u s 

m a -1 

<5z s 

m 

0.50 a s 

6 u s 
m a -1 

SZs 

m 

0.75 a s 

6 u s 
m a' 1 

2 

32 (5%) 

-11 (-9%) 

2 (0%) 

-27 (-24%) 

-101 (-15%) 

-42 (-36%) 

3 

48 (7%) 

1 7 (1 5%) 

30 (4%) 

4 (3%) 

-11 (-2%) 

-1 9 (-1 7%) 

4 

59 (9%) 

40 (34%) 

45 (7%) 

28 (24%) 

25 (4%) 

12 (11%) 
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Fig. 7. Modeled time-space distribution of basal sliding velocity, u b, 
along the terminal 60km of the Sermeq Avannarleq flowline. 
Vertical dashed lines indicate the locations of JAR2, JAR1 and Swiss 
Camp (SC). Color bar saturates at lOOmaW 


satisfactory annual basal sliding cycle (Fig. 7). We interpret 
this as suggesting that the dhjdt term contains important 
information that modulates basal sliding. This suggests that 
sliding rules that are dependent on basal stress and absolute 
head may overlook important variables that are correlated 
with dh e /dt, particularly those related to transient subglacial 
transmissivity or higher-order features of subglacial hydro- 
logic geometry (e.g. the r? c (x) 'conduit spacing' parameter 
used in this study). Recent observations of hysteresis in the 
ratio of inferred subglacial water storage to sliding velocity 
over the course of a melt season support this notion (Howat 
and others, 2008a). Thus, we suggest that an improvement 
towards achieving a physically based sliding rule would be 
to combine an absolute head/basal stress rule that modulates 
background (i.e. winter) sliding velocities (which are 
prescribed in this study) with a dh e /dt-type parameterization 
that honors the empirical subtlety that sliding reaches a 
maximum when dh e /dt reaches a maximum (rather than 
when h e reaches a maximum). 

4.1 . Modeled speed-up and slowdown events 

There are discrepancies in the absolute magnitude and 
timing between modeled and observed velocity maxima and 
minima at all three stations (Fig. 8). As our interest lies in 
producing an annual velocity cycle using an annual 
hydrologic cycle, we focus the following discussion primar- 
ily on temporal discrepancies. We speculate that the 
relatively damped annual velocity cycle observed at JAR1 
is due to a relatively dampened annual hydrology cycle. 
While the hydrology model forces subglacial water to flow 
over the bedrock high at km 30, in reality subglacial water 
likely flows around this bedrock high in the /direction. Thus, 
the hydrology model likely overestimates 9h e /c>fvalues along 
the km 25-40 portion of the flowline, for which the scaling 
parameter k(x) compensates. In regard to the absolute 
magnitude of the annual cycle, we note that at any location 
along the flowline an infinite combination of the parameters 
k(x) and m may be used to scale dhjdtto a desired value. We 
also note that the modeled basal sliding velocity cycles 
exhibit abrupt transitions between background sliding and 
speed-up and slowdown events. These abrupt transitions are 
due to the imposition of a relatively high threshold of rate of 
change of glacier water storage to transition between the 



Fig. 8. Modeled and observed basal sliding velocity versus day of 
year at Swiss Camp (SC), JAR1 and JAR2. Discrepancies in the 
timing of the summer speed-up and fall slowdown events are 
denoted with dashed lines. 


three sliding phases {\dh e /dt\ >0.25 md" 1 ). In reality, this 
threshold likely varies in time and space. 

We take uncertainty in the timing discrepancy as 
equivalent to the temporal resolution of the velocity data 
(i.e. ±10 days). Modeled maximum velocity precedes 
observed maximum velocity by ~20±10 days at Swiss 
Camp and ~1 8 ± 1 0 days at JAR1 (Fig. 8). It is more difficult 
to assess the timing discrepancy of the summer speed-up 
event at JAR2, where an apparently premature onset of the 
modeled slowdown event truncates the speed-up event, but 
the discrepancy appears to be <10± 10 days. The relatively 
early onset of the modeled speed-up event can be attributed 
to a combination of: (1) incorrect timing of the prescribed 
meltwater input and (2) the assumption embedded in the 
hydrology model that meltwater produced at the surface is 
immediately routed to the glacier interior, where it raises the 
englacial water table. This second assumption assumes that 
there is no temporary supraglacial meltwater storage (e.g. 
within a saturated snowpack or in ponded water) and that 
supra- and englacial travel times are negligible. In reality, 
temporary supraglacial meltwater storage and travel time 
can delay the initial meltwater pulse from reaching the 
englacial water table for several weeks after the onset of 
melt (Fountain and Walder, 1 998; Flowers and Clarke, 2002; 
Jansson and others, 2003). This lag would be expected to be 
greater near the equilibrium line (i.e. at Swiss Camp) than at 
the terminus (i.e. at JAR2). In addition to temporary 
supraglacial meltwater storage, englacial transfer time can 
be expected to range over two orders of magnitude as a 
consequence of the morphology of the englacial hydrologic 
system (i.e. moulin- versus crevasse-type drainage; Colgan 
and others, 2011b; McGrath and others, 2011). A more 
detailed treatment of supra- and englacial meltwater routing 
in the hydrology model would likely reduce the discrep- 
ancies between modeled and observed velocity maxima. 

The discrepancy between modeled and observed velocity 
minima is more variable; the model produces a significantly 
earlier onset at JAR2 (~65±10 days), an extreme delay at 
Swiss Camp (~110±10 days) and a reasonable match at 
JAR1 (1 5 ± 1 0 days; Fig. 8). We believe that this reflects the 
difficulty in capturing the behavior of subglacial conduits in 
an along-flowline 1 -D hydrology model. The timing of the fall 
slowdown event depends on the timing of the negative (or 
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Distance upstream (km) 

Fig. 9. Observed day of year of maximum (/ max ; red) and minimum 
(y min ; blue) ice surface velocity at the GC-Net stations versus 
distance upstream. Vertical whiskers denote ±10 days uncertainty 
in /max and / min at each station. Solid lines denote least-squares 
linear fit. Dashed lines denote 95% confidence bounds for the 
linear fit. The grey shading denotes the 95% confidence envelope 
for the convergence of y max and y min . The black dot denotes day 223 
and upstream distance km 69.5. 


decreasing) glacier water storage ( dhjdt ) phase. This phase is 
initiated by the upstream propagation of a nickpoint in 
englacial water table elevation that corresponds to the 
opening of efficient subglacial conduits that in turn lower 
glacier water storage (Animation 1; e.g. Kessler and 
Anderson, 2004). A two-dimensional (2-D) (xy) hydrological 
model, which allows water to flow both parallel and perpen- 
dicular to the ice dynamic flowline, would be inherently 
more realistic in propagating changes in englacial water table 
elevation upstream from various outlets at the ice margin. 
Because the timing of the slowdown event is controlled by 
the evolution of conduit sizes, the fall slowdown event is 
more dependent on an accurate representation of subglacial 
conduit dynamics than is the spring speed-up event. The 
spring speed-up event initiates as long as conduits have 
collapsed to their minimum radii (i.e. 'closed'), which is a 
relatively simple geometry to capture given the lengthy 
winter period over which this geometry is achieved. 

4.2. Upstream limit of seasonal basal sliding 

The observed dates of summer maximum (y max ) and fall 
minimum (y min ) velocity at JAR2, JAR1 and Swiss Camp 
suggest that the timing of the fall slowdown event is more 
synchronous than the timing of the summer speed-up event 
along the flowline (Fig. 9). GPS observations indicate that 
the summer speed-up event 'propagates' (or is transmitted) 
upstream at 1.2kmd _1 (with a 95% confidence interval of 
0.8-2. 6 km d -1 ), while the fall slowdown event propagates 
downstream at 1.7kmd _1 (with a 95% confidence interval 
of 1.0-5. 2 kmd -1 ). We interpret this as suggesting that the 
timing of the speed-up event depends on a meteorological 
forcing (i.e. the onset of melt or a critical surface ablation 
threshold), while the timing of the slowdown event is 
dependent on the development of efficient subglacial 
transmission capacity (which is more likely to be synchro- 
nous over the flowline). We speculate that the theoretical 
upstream limit to which the annual basal sliding cycle can 
propagate (i.e. the distance inland of which the ice sheet 
does not experience an annual velocity cycle) is defined by 
the upstream convergence of the dates of maximum and 
minimum velocity. The 2005 and 2006 GPS data suggest 



Distance upstream (km) 

Fig. 10. Total driving stress, r, and depth-averaged longitudinal 
coupling stress, t' xx , along the terminal 50km of the flowline at 
dynamic equilibrium (i.e. post-spin-up). Dashed lines represent 
±10% of total driving stress. 


that this occurs at ~69.5 km upstream on day of year ~223. 
This distance upstream corresponds to the ~1350m ice 
elevation contour, which is slightly above the regional ELA 
over the period 1996-2006 (~1250m; Fausto and others, 
2009). The lower limit of this convergence position, 
delineated by the 95% confidence envelope, is ~55 km 
upstream. GPS velocity observations at Up50 (km 98.5 at 
69.75° N, 48.14° W in 2008), which show no evidence of an 
annual velocity cycle (Colgan and others, 2009), provide an 
upstream limit for the annual velocity cycle. 

4.3. Importance of longitudinal coupling 

Although in situ GPS data indicate that Swiss Camp exhibits 
an annual velocity cycle, these data cannot indicate whether 
this is due to: (1) local meltwater-induced acceleration 
(Zwally and others, 2002; Bartholomew and others, 201 0) or 
(2) an annual velocity cycle originating downstream of Swiss 
Camp that is propagated upstream via longitudinal coupling. 
This downstream annual velocity cycle may be due to either 
lower-elevation meltwater-induced acceleration (Price and 
others, 2008) or the annual tidewater calving cycle (Howat 
and others, 2008b; Joughin and others, 2008). A previous 
study has suggested that a 1 0-20% velocity increase at Swiss 
Camp can be achieved by a roughly 100% velocity increase 
initiated ~12km downstream from Swiss Camp that is 
propagated upstream through longitudinal coupling (Price 
and others, 2008). This previous study used a 2-D (cross- 
sectional) ice-flow model in which basal sliding was par- 
ameterized to occur through deformation of a fluid layer 
several meters thick underlying the ice. This basal fluid layer 
was assumed to have an effective viscosity of ~2.8 x 1 0 4 Pa a. 
For comparison, typical flow-law parameter values 
(i.e. A« 10 -16 Pa -3 a -1 ) represent an effective viscosity of 
~2.5 x 1 0 6 Pa a when H= 500 m, dzjdx= 0.01 and n = 3. In 
addition to a vertical structure that is conducive to transmit- 
ting longitudinal coupling, the magnitude of the assumed 
downstream seasonal velocity perturbation (i.e. 100% or 
doubling) greatly exceeds that recorded by JAR1 station 
1 3.5 km downstream of Swiss Camp. 

In the present study, the 1-D (depth-integrated) ice-flow 
model suggests that the absolute longitudinal coupling stress, 
|r' vx |; Eqn (6), is only significant (defined here as >10% of 
total driving stress) along the terminal ~6 km of the flowline 
(Fig. 10). Upstream of the icefall at ~6 km, where observed 
ice thickness decreases to <300 m, the absolute longitudinal 
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Model time: 5.996 



Distance upstream (km) 

Animation 2. Animation of the terminus perturbation simulation, (a) Total driving stress, r, and depth-averaged longitudinal coupling stress, 
r' xx . (b) Bedrock elevation (brown), dynamic equilibrium ice geometry at year 0 (grey line) and transient ice geometry (black line), with 
contour shading to denote ice surface velocity, u s (color bar saturates at 200 ma -1 ). The vertical dashed red line denotes the location of the 
first-type boundary condition imposed to represent a catastrophic terminus perturbation, (c) Prescribed background basal sliding velocity. 
Vertical dashed lines denote the locations of JAR2, JAR1 and Swiss Camp (SC). Model time (in years) given as relative to the end of spin-up. 
Full movie available at htt p://www.igsoc.org/hyp er link/11l081 Animation2.mov 


coupling stresses seldom exceed 10% of total driving stress. 
In addition, any perturbation to the tidewater tongue would 
also experience rapid radial diffusion (in the x-y plane) with 
distance inland. This inference of minimal inland coupling 
stresses fits the theoretical notion that coupling stresses are 
typically important only in the terminal few km of ice-sheet 
flowlines (where ice thickness becomes small). In these 
terminal regions the magnitude of coupling stresses can be 
equivalent to, or exceed, driving stresses (Van der Veen, 
1 987). Where longitudinal coupling stresses are insignificant 
along the Sermeq Avannarleq flowline, the forces governing 
ice flow can be assumed to be local in nature. 

The inability of a seasonal terminus perturbation (cf. 
Howat and others, 2008b; Joughin and others, 2008) to 
propagate upstream to produce the ~55% seasonal velocity 
acceleration observed at Swiss Camp (from u w = 1 1 3 m a -1 to 
u m ax= 1 75 m a -1 ) can also be demonstrated by a simple 
numerical simulation. In this fully transient simulation, we 
spin-up the ice-flow model for 1000 years under the £=3 
and 50% of contemporary a s scenario. To isolate potential 
changes in ice dynamics, we disable hydrological coupling 
and prescribe temporally constant background basal sliding 
along the flowline. At 1 year after spin-up, we impose a first- 
type (specified head) Dirichlet boundary condition of 
H=0m at km 8. This boundary condition, which instantly 
removes the terminal several km of the flowline, is meant to 
represent a catastrophic terminus perturbation. Despite 
imposing an unprecedented terminus perturbation that is 
an order of magnitude larger than the annual advance/retreat 
cycle, longitudinal coupling stresses inland of JAR2 (km 14) 
fail to exhibit any significant change in the 4 years following 
the perturbation (Animation 2). The ice acceleration is 
propagated up to Swiss Camp (km 46) not by longitudinal 
stresses but by changes in ice geometry (i.e. surface slope 
steepening) over ~2.75 years (Fig. 11). The modeled 


upstream propagation rate (~17kma -1 ) is over an order of 
magnitude slower than the upstream propagation rate of the 
summer speed-up event inferred from the GPS observations 
(~1.2kmd -1 ). In addition, upon reaching Swiss Camp the 
modeled velocity perturbation is only ~5% of the modeled 
dynamic equilibrium winter velocity at Swiss Camp, which 
is an order of magnitude less than the observed ~55% 
seasonal acceleration. Thus, it is unlikely that terminus 
perturbations associated with the annual retreat/advance 
cycle (Howat and others, 2008b; Joughin and others, 2008) 
are propagated upstream by longitudinal coupling (Price and 
others, 2008) to produce the seasonal velocity variations 
observed inland of JAR2 (km 14). Instead, we interpret the 
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Fig. 11. Time-space distribution of ice surface velocity, u s , in the 
terminus perturbation simulation (color bar saturates at 200 ma -1 ). 
Labeled black contours represent the relative magnitude of the 
velocity perturbation (relative to dynamic equilibrium velocity in 
year 0) resulting from the catastrophic terminus perturbation event at 
year 1 . Vertical dashed lines denote the locations of JAR2, JAR1 and 
Swiss Camp (SC). Model time given as relative to the end of spin-up. 
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annual velocity cycle at upstream sites to reflect local 
mismatches between glacier water inputs and outputs (i.e. 
local melt-induced acceleration; Zwally and others, 2002; 
Kessler and Anderson, 2004). 

5. SUMMARY REMARKS 

We examined the annual glaciohydrology cycle in the 
ablation zone of the Sermeq Avannarleq flowline. We 
coupled a 1-D (depth-integrated) ice-flow model to a 
previously published 1-D (depth-integrated) hydrology 
model via a semi-empirical and site-specific basal sliding 
rule. This sliding rule imposes seasonal perturbations to the 
background basal sliding velocity that are dependent on rate 
of change of glacier water storage. Following a lOOOyear 
spin-up, the ice-flow model produces dynamic equilibrium 
ice geometry and velocity fields that compare well with 
observations. After spin-up, the coupled model reproduces 
the broad features observed in the annual basal sliding cycle 
in the terminal 60 km of the flowline, namely: (1 ) prescribed 
background basal sliding during the winter; (2) a summer 
speed-up event; and (3) a fall slowdown event and return to 
winter velocities. While we have put forth a plausible rule 
that connects sliding velocity to the state of the glacier 
hydrologic system, there remains a significant challenge to 
develop a sliding rule that is more firmly based upon first 
principles (i.e. no free parameters and based on physical 
mechanisms by which basal water produces sliding). This 
requires both a proper characterization of all components of 
the glacier hydrologic system, from the spatial-temporal 
pattern of melt generation to the complex evolution of en- 
and subglacial transport and storage, and a physics-based 
connection between the state of the hydrologic system and 
the basal sliding. 

GPS observations of ice surface velocity at three stations 
during 2005 and 2006 suggest the annual velocity cycle 
propagates as far upstream as ~70km. We examined the 
relative magnitude of driving and coupling stresses and 
performed a simple simulation to assess the possible contri- 
bution of an upstream propagation of a terminus perturbation 
through longitudinal coupling to the annual velocity cycle. 
We find that an extreme terminus perturbation is unlikely to 
influence velocities upstream of JAR2 (km 14) on the seasonal 
timescale. Thus, we suggest the annual ice velocity cycle 
along the majority of the flowline (km 14-70) is instead 
attributable to the evolution of the glaciohydrologic system in 
response to meltwater inputs. Following previous alpine 
studies, we suggest this local acceleration is due to variations 
in basal sliding that are governed by variations in the rate of 
change of glacier water storage, due to local mismatches 
between surface meltwater input and the ability of the 
subglacial hydrologic system to transmit water. 
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APPENDIX 

Variable notation 

7 Ablation gradient 

Pi Density of ice (kg rrT 3 ) 

t Total driving stress (Pa) 

t' xx Longitudinal stress (Pa) 

A Flow-law parameter (Pa -3 a -1 ) 

E Wisconsin enhancement factor 

F Englacial hydrology system entry fraction 
H Ice thickness (m) 

K Sliding rule coefficient (m 1 / 3 d 2/3 ) 

N Glen law exponent 

P\ Ice pressure (Pa) 

P w Water pressure (Pa) 

Q Ice discharge per unit width (m 2 a _1 ) 

Q e Creep activation energy of ice (kj mol ') 

R Ideal gas constant (J mob 1 K” 1 ) 

T Ice temperature (K) 

ab Annual basal ablation (m) 

a s Annual surface ablation (m) 

a s Surface ablation rate (ma -1 ) 

b Annual mass balance (m) 
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Cb Annual basal accumulation (m) 

c s Annual surface accumulation (m) 

4nax Duration of summer speed-up event (days) 
c( min Duration of fall slowdown event (days) 
g Gravitational acceleration (m s -2 ) 

h e Englacial water table elevation (or head) (m) 

/ Day of year (days) 

/max Day of year of summer maximum velocity (days) 

/min Day of year of fall minimum velocity (days) 
m Sliding rule exponent 

n c Subglacial conduits per meter (m -1 ) 
t Given time (years) 

Ub Basal sliding velocity (ma' 1 ) 


Ubo Background basal sliding velocity (ma -1 ) 

Ud Deformational velocity (m a -1 ) 

Umax Ice surface summer maximum velocity (ma -1 ) 

Umin Ice surface fall minimum velocity (ma -1 ) 
u s Ice surface velocity (m a -1 ) 

u w Ice surface winter velocity (ma -1 ) 
x Given distance upstream (m) 

z Given elevation (m) 

Zb Bedrock elevation (m) 

z s Ice surface elevation (m) 

All heads and elevations are in reference to sea level as 
datum. 
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